# # # # # # # # # # # # # # # # # # #
# Autor: Tiago Chedraoui Silva
# RA: 082941
# Descricao: http://tinyurl.com/3qa2yze
# # # # # # # # # # # # # # # # # # # 

init <- function(){
  # salva linha em vetores
  data <- array((scan("digitos-sem-classe.txt", what = integer(0))),, dim=c(64,3823))

  # Inicializacao dos 10 centros 
  centers <-array(dim=c(64,10))
  centers[,1] <-  data[,1]  # Numero 0
  centers[,2] <-  data[,12] # Numero 1
  centers[,3] <-  data[,26] # Numero 2
  centers[,4] <-  data[,15] # Numero 3
  centers[,5] <-  data[,4]  # Numero 4
  centers[,6] <-  data[,7]  # Numero 5
  centers[,7] <-  data[,59]  # Numero 6
  centers[,8] <-  data[,3]  # Numero 7
  centers[,9] <-  data[,10] # Numero 8
  centers[,10] <- data[,24] # Numero 9

  centers<- myKmean(data,centers)
  drawCenters <- CenterToInt(centers)
  printCenters(drawCenters)
  # Verificacao do metodo:imprime centros em figura
  CentersToPGM(drawCenters,"img/ponto")

# Inicializacao dos 10 centros com outros pontos
  centers <-array(dim=c(64,10))
  centers[,1] <-  data[,2]  # Numero 0
  centers[,2] <-  data[,70] # Numero 1
  centers[,3] <-  data[,3782] # Numero 2
  centers[,4] <-  data[,75] # Numero 3
  centers[,5] <-  data[,95]  # Numero 4
  centers[,6] <-  data[,107]  # Numero 5
  centers[,7] <-  data[,3655]  # Numero 6
  centers[,8] <-  data[,3741]  # Numero 7
  centers[,9] <-  data[,89] # Numero 8
  centers[,10] <- data[,41] # Numero 9

  centers<- myKmean(data,centers)
  drawCenters <- CenterToInt(centers)
  printCenters(drawCenters)
  #Verificacao do metodo:imprime centros em figura
  CentersToPGM(drawCenters,"img/precisos")
  

  # Melhorando o valor do centros iniciais
  centers<-myKmean(data,centers)
  centers <- CenterToInt(centers)
  printCenters(centers)

  # Verificacao do metodo:imprime centros em figura
  CentersToPGM(centers,"img/ponto_melhor")
  
  # Verificacao do metodo:Semente randomica
  for(aux in 1:10){
    centers[,aux] <- c(sample(0:16, 64, replace =TRUE))
  }
  centers<-myKmean(data,centers)
  printCenters(centers)

  centers <- CenterToInt(centers)
  CentersToPGM(centers,"img/rnd")

  # INICIO PARTE 5
  # Dividimos os pontos em grupos de 10  Transformamos o centroide de cada grupo na origem
  # Calculosmo a matrix de covariancia S e calculamos o
  # dois componentes principais(PCA)
  # Saida para cada grupo: m +- s1*v1 e m +- s2*v2
  r <- array(dim=c(10,3823))
  for (i in 1:3823){
    r[,i]= findGroup(data[,i],centers)
  }
  
  x<-{}
  y<-{}
  for(k in 1:10){
    write(paste("== GRUPO",k,"===="),"pca.txt",append = TRUE)
    grupok<-groupDiv(data,k,r)
    grupok<-centroOrigem(grupok,centers[,k])
    S<-array()
    S<-CoVarMatrix(grupok)
    v<-PrincipalComp(S,centers[,k],grupok) 
    x<-c(x,v[1])
    y<-c(y,v[1])
  }
  plot(x,y)
  # FIM PARTE 5
}
# Fim Parte 1 - Inicio Parte 2

myKmean <- function(x,centers){
  count=0
  while(1){
    # Parte um do algoritmo: Definicao dos grupos
    # Cria um vetor r para cada ponto
    # r[k,i] = 1 se i pertence a k, 0 caso contrario
    r <- array(dim=c(10,3823))
    for (i in 1:3823){
      r[,i]= findGroup(x[,i],centers)
    }
    # Parte dois do algoritmo: Atualizacao de centros

    # Sera que nosso algortimo convergiu?
    # Guarde-se centros atuais para posterior comparacao
    oldCenters = centers
    # para cada centro calcula numero de pontos no grupo
    pointsInGroup <- array(data=0,dim=c(10))
    # para cada centro calcula soma dos valores dos pontos no grupo
    sumInGroup <-array(data=0,dim=c(64,10))

    #Calculo de novo centro
    for (k in 1:10){
      for (i in 1:3823){
        pointsInGroup[k]=pointsInGroup[k]+(r[k,i])
        sumInGroup[,k]=sumInGroup[,k]+(r[k,i])*x[,i]
      }
      if(pointsInGroup[k]>0){
        centers[,k]=sumInGroup[,k]/pointsInGroup[k]
      }
    }

    #Fim da convergencia?
    if(sum(oldCenters-centers)==0){
      print ("ACHEI OS CENTROS")
      print(count)
      break 
    }
    else if(count>100){
      break
    }
    count <- count + 1
    drawCenters <- CenterToInt(centers)
    printCenters(drawCenters)
    #Verificacao do metodo:imprime centros em figura
    CentersToPGM(drawCenters,paste("img/",count))

    }
  return (centers)
}

# Parte 5
# Dado o conjunto de dados devo altera a origem deles
# novo x = x - c
centroOrigem<-function(x,centers){
  
  # Parte um do algoritmo: Definicao dos grupos
  
  for (i in 1:(length(x)/64)){
    x[,i]=x[,i]-centers
  }
    return (x)
}

groupDiv<- function(x,number,r){
  group <-array(dim=c(64,1))
  for (i in 1:3823){
    if(r[number,i]==1){
      group<-cbind(group,x[,i])
    }
  }
  return (group[,-1])
}


PrincipalComp<-function(S,center,grupo){
  dc <- svd(S)
  #Vetor diagonal da matriz
  D<-dc$d
  # Matrizes da decomposicao
  U<-dc$u
  V<-dc$v
  # primeiros dois principais componentes
  u1<-dc$u[,1]
  u2<-dc$u[,2]

  u1Mean<-sum(u1)/length(u1)
  u1Sigma<- ((u1[1]-(u1Mean))^2)
  for(i in 2:(64)){
    u1Sigma<- u1Sigma+((u1[i]-(u1Mean))^2)
  }
  u1Sigma<-sqrt(u1Sigma/63)

  u2Mean<-sum(u2)/length(u2)
  u2Sigma<- ((u2[1]-(u2Mean))^2)
  for(i in 2:(64)){
    u2Sigma<- u2Sigma+((u2[i]-(u2Mean))^2)
  }
  u2Sigma<-sqrt(u2Sigma/63)

  x= sum((center+u1Sigma*u1)*u1)
  y= sum((center+u2Sigma*u2)*u2)
  print(x)
  print(y)
  write("m + s1*v1", "pca.txt", sep = "\t",append = TRUE)
  write(center+u1Sigma*u1, "pca.txt", sep = "\t",append = TRUE)
  write("m - s1*v1", "pca.txt", sep = "\t",append = TRUE)
  write(center-u1Sigma*u1, "pca.txt", sep = "\t",append = TRUE)
  write("m + s2*v2", "pca.txt", sep = "\t",append = TRUE)
  write(center+u2Sigma*u2, "pca.txt", sep = "\t",append = TRUE)
  write("m - s2*v2", "pca.txt", sep = "\t",append = TRUE)
  write(center-u2Sigma*u2, "pca.txt", sep = "\t",append = TRUE)
  v<- c(x,y)
  return (v)
}
  
CoVarMatrix <- function(x){
  S <-array(data=0,dim=c(64,64))
  # para cada centro calcula numero de pontos no grupo
  # vou encontrar os pontos no grupo e aplicar
  S = S + x%*%t(x)
  return (S)
}

#Create PGM from a matrix
# Entrada vetor
# Saida: Arquivo PGM 
MatrixToPGM <- function(vect){
  write("P2", "data.pgm",append = FALSE)
  write("8 8", "data.pgm",append = TRUE)
  write("16", "data.pgm",append = TRUE)
  write(vect, "data.pgm", sep = "\t",append = TRUE)
}

# Create PGM from a matrix
# Entrada vetor
# Saida: Arquivo PGM
CentersToPGM <- function(vect,name){
  for (i in 1:10){
    write("P2", paste(name,"center_",i,".pgm"),append = FALSE)
    write("8 8", paste(name,"center_",i,".pgm"),append = TRUE)
    write("16", paste(name,"center_",i,".pgm"),append = TRUE)
    write(vect[,i], paste(name,"center_",i,".pgm"), sep = "\t",append = TRUE)
  }
}

# Imprime centroide final
printCenters <- function(centers){
  for(aux in 1:10){
    print(paste("== Centroide do numero:",aux,"=="))
    print (centers[,aux])
  }
}

# Calcula distancia eucliadiana entre dois pontos em dimensao n
# Entrada: Dois vetores X e Y de dimensao n
# Saida distancia euclidiana
euclDist <- function(p1,p2){
    return(sqrt(sum((p1 -p2)^2)))
}

# Procura em qual grupo o ponto deve estar
# Entrada: ponto 
# Saida vetor r que vale:
# 1 para o grupo ao qual o ponto pertence 
# 0 para os outros
findGroup <- function(p,centers){
  min = Inf
  min_idx=1
  for (i in 1:10){
    dist=euclDist(p,centers[,i])
  
    if(dist<min){
      min=dist
      min_idx=i
    }
  }
  
  r<- c(seq(0,0,length=10))
  r[min_idx]=1
  return (r)
}

# transforma numeros fracionarios em inteiros
CenterToInt<- function(centers){
  for(aux in 1:10){
    for(i in 1:64){
      centers[i,aux] <- as.integer(centers[i,aux])
    }
  }
  return (centers)
}

CentersToPGM <- function(vect,name){
  for (i in 1:3823){
    write("P2", paste(name,i,".pgm"),append = FALSE)
    write("8 8", paste(name,i,".pgm"),append = TRUE)
    write("16", paste(name,i,".pgm"),append = TRUE)
    write(vect[,i], paste(name,i,".pgm"), sep = "\t",append = TRUE)
  }
}

